home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qmomof.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-04-17  |  8.4 KB  |  390 lines
  1. /* integration/qmomof.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_integration.h>
  23. #include <gsl/gsl_errno.h>
  24.  
  25. static void
  26. compute_moments (double par, double * cheb);
  27.  
  28. static int
  29. dgtsl (size_t n, double *c, double *d, double *e, double *b);
  30.  
  31. gsl_integration_qawo_table *
  32. gsl_integration_qawo_table_alloc (double omega, double L, 
  33.                   enum gsl_integration_qawo_enum sine,
  34.                   size_t n)
  35. {
  36.   gsl_integration_qawo_table *t;
  37.   double * chebmo;
  38.  
  39.   if (n == 0)
  40.     {
  41.       GSL_ERROR_VAL ("table length n must be positive integer",
  42.             GSL_EDOM, 0);
  43.     }
  44.  
  45.   t = (gsl_integration_qawo_table *)
  46.     malloc (sizeof (gsl_integration_qawo_table));
  47.  
  48.   if (t == 0)
  49.     {
  50.       GSL_ERROR_VAL ("failed to allocate space for qawo_table struct",
  51.             GSL_ENOMEM, 0);
  52.     }
  53.  
  54.   chebmo = (double *)  malloc (25 * n * sizeof (double));
  55.  
  56.   if (chebmo == 0)
  57.     {
  58.       free (t);
  59.       GSL_ERROR_VAL ("failed to allocate space for chebmo block",
  60.             GSL_ENOMEM, 0);
  61.     }
  62.  
  63.   t->n = n;
  64.   t->sine = sine;
  65.   t->omega = omega;
  66.   t->L = L;
  67.   t->par = 0.5 * omega * L;
  68.   t->chebmo = chebmo;
  69.  
  70.   /* precompute the moments */
  71.  
  72.   { 
  73.     size_t i;
  74.     double scale = 1.0;
  75.  
  76.     for (i = 0 ; i < t->n; i++)
  77.       {
  78.         compute_moments (t->par * scale, t->chebmo + 25*i);
  79.         scale *= 0.5;
  80.       }
  81.   }
  82.  
  83.   return t;
  84. }
  85.  
  86. int
  87. gsl_integration_qawo_table_set (gsl_integration_qawo_table * t,
  88.                     double omega, double L,
  89.                     enum gsl_integration_qawo_enum sine)
  90. {
  91.   t->omega = omega;
  92.   t->sine = sine;
  93.   t->L = L;
  94.   t->par = 0.5 * omega * L;
  95.  
  96.   /* recompute the moments */
  97.  
  98.   { 
  99.     size_t i;
  100.     double scale = 1.0;
  101.  
  102.     for (i = 0 ; i < t->n; i++)
  103.       {
  104.         compute_moments (t->par * scale, t->chebmo + 25*i);
  105.         scale *= 0.5;
  106.       }
  107.   }
  108.  
  109.   return GSL_SUCCESS;
  110. }
  111.  
  112.  
  113. int
  114. gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * t,
  115.                        double L)
  116. {
  117.   /* return immediately if the length is the same as the old length */
  118.  
  119.   if (L == t->L)
  120.     return GSL_SUCCESS;
  121.  
  122.   /* otherwise reset the table and compute the new parameters */
  123.  
  124.   t->L = L;
  125.   t->par = 0.5 * t->omega * L;
  126.  
  127.   /* recompute the moments */
  128.  
  129.   { 
  130.     size_t i;
  131.     double scale = 1.0;
  132.  
  133.     for (i = 0 ; i < t->n; i++)
  134.       {
  135.         compute_moments (t->par * scale, t->chebmo + 25*i);
  136.         scale *= 0.5;
  137.       }
  138.   }
  139.  
  140.   return GSL_SUCCESS;
  141. }
  142.  
  143.  
  144. void
  145. gsl_integration_qawo_table_free (gsl_integration_qawo_table * t)
  146. {
  147.   free (t->chebmo);
  148.   free (t);
  149. }
  150.  
  151. static void
  152. compute_moments (double par, double *chebmo)
  153. {
  154.   double v[28], d[25], d1[25], d2[25];
  155.  
  156.   const size_t noeq = 25;
  157.   
  158.   const double par2 = par * par;
  159.   const double par4 = par2 * par2;
  160.   const double par22 = par2 + 2.0;
  161.  
  162.   const double sinpar = sin (par);
  163.   const double cospar = cos (par);
  164.  
  165.   size_t i;
  166.  
  167.   /* compute the chebyschev moments with respect to cosine */
  168.  
  169.   double ac = 8 * cospar;
  170.   double as = 24 * par * sinpar;
  171.  
  172.   v[0] = 2 * sinpar / par;
  173.   v[1] = (8 * cospar + (2 * par2 - 8) * sinpar / par) / par2;
  174.   v[2] = (32 * (par2 - 12) * cospar
  175.       + (2 * ((par2 - 80) * par2 + 192) * sinpar) / par) / par4;
  176.  
  177.   if (fabs (par) <= 24)
  178.     {
  179.       /* compute the moments as the solution of a boundary value
  180.          problem using the asyptotic expansion as an endpoint */
  181.       
  182.       double an2, ass, asap;
  183.       double an = 6;
  184.       size_t k;
  185.  
  186.       for (k = 0; k < noeq - 1; k++)
  187.     {
  188.       an2 = an * an;
  189.       d[k] = -2 * (an2 - 4) * (par22 - 2 * an2);
  190.       d2[k] = (an - 1) * (an - 2) * par2;
  191.       d1[k + 1] = (an + 3) * (an + 4) * par2;
  192.       v[k + 3] = as - (an2 - 4) * ac;
  193.       an = an + 2.0;
  194.     }
  195.  
  196.       an2 = an * an;
  197.  
  198.       d[noeq - 1] = -2 * (an2 - 4) * (par22 - 2 * an2);
  199.       v[noeq + 2] = as - (an2 - 4) * ac;
  200.       v[3] = v[3] - 56 * par2 * v[2];
  201.  
  202.       ass = par * sinpar;
  203.       asap = (((((210 * par2 - 1) * cospar - (105 * par2 - 63) * ass) / an2
  204.         - (1 - 15 * par2) * cospar + 15 * ass) / an2 
  205.            - cospar + 3 * ass) / an2 
  206.           - cospar) / an2;
  207.       v[noeq + 2] = v[noeq + 2] - 2 * asap * par2 * (an - 1) * (an - 2);
  208.  
  209.       dgtsl (noeq, d1, d, d2, v + 3);
  210.  
  211.     }
  212.   else
  213.     {
  214.       /* compute the moments by forward recursion */
  215.       size_t k;
  216.       double an = 4;
  217.  
  218.       for (k = 3; k < 13; k++)
  219.     {
  220.       double an2 = an * an;
  221.       v[k] = ((an2 - 4) * (2 * (par22 - 2 * an2) * v[k - 1] - ac)
  222.           + as - par2 * (an + 1) * (an + 2) * v[k - 2]) 
  223.         / (par2 * (an - 1) * (an - 2));
  224.       an = an + 2.0;
  225.     }
  226.     }
  227.  
  228.  
  229.   for (i = 0; i < 13; i++)
  230.     {
  231.       chebmo[2 * i] = v[i];
  232.     }
  233.  
  234.   /* compute the chebyschev moments with respect to sine */
  235.  
  236.   v[0] = 2 * (sinpar - par * cospar) / par2;
  237.   v[1] = (18 - 48 / par2) * sinpar / par2 + (-2 + 48 / par2) * cospar / par;
  238.  
  239.   ac = -24 * par * cospar;
  240.   as = -8 * sinpar;
  241.  
  242.   if (fabs (par) <= 24)
  243.     {
  244.       /* compute the moments as the solution of a boundary value
  245.          problem using the asyptotic expansion as an endpoint */
  246.  
  247.       size_t k;
  248.       double an2, ass, asap;
  249.       double an = 5;
  250.  
  251.       for (k = 0; k < noeq - 1; k++)
  252.     {
  253.       an2 = an * an;
  254.       d[k] = -2 * (an2 - 4) * (par22 - 2 * an2);
  255.       d2[k] = (an - 1) * (an - 2) * par2;
  256.       d1[k + 1] = (an + 3) * (an + 4) * par2;
  257.       v[k + 2] = ac + (an2 - 4) * as;
  258.       an = an + 2.0;
  259.     }
  260.       
  261.       an2 = an * an;
  262.  
  263.       d[noeq - 1] = -2 * (an2 - 4) * (par22 - 2 * an2);
  264.       v[noeq + 1] = ac + (an2 - 4) * as;
  265.       v[2] = v[2] - 42 * par2 * v[1];
  266.  
  267.       ass = par * cospar;
  268.       asap = (((((105 * par2 - 63) * ass - (210 * par2 - 1) * sinpar) / an2
  269.         + (15 * par2 - 1) * sinpar
  270.         - 15 * ass) / an2 - sinpar - 3 * ass) / an2 - sinpar) / an2;
  271.       v[noeq + 1] = v[noeq + 1] - 2 * asap * par2 * (an - 1) * (an - 2);
  272.  
  273.       dgtsl (noeq, d1, d, d2, v + 2);
  274.  
  275.     }
  276.   else
  277.     {
  278.       /* compute the moments by forward recursion */
  279.       size_t k;
  280.       double an = 3;
  281.       for (k = 2; k < 12; k++)
  282.     {
  283.       double an2 = an * an;
  284.       v[k] = ((an2 - 4) * (2 * (par22 - 2 * an2) * v[k - 1] + as)
  285.           + ac - par2 * (an + 1) * (an + 2) * v[k - 2]) 
  286.         / (par2 * (an - 1) * (an - 2));
  287.       an = an + 2.0;
  288.     }
  289.     }
  290.  
  291.   for (i = 0; i < 12; i++)
  292.     {
  293.       chebmo[2 * i + 1] = v[i];
  294.     }
  295.  
  296. }
  297.  
  298. static int
  299. dgtsl (size_t n, double *c, double *d, double *e, double *b)
  300. {
  301.   /* solves a tridiagonal matrix A x = b 
  302.      
  303.      c[1 .. n - 1]   subdiagonal of the matrix A
  304.      d[0 .. n - 1]   diagonal of the matrix A
  305.      e[0 .. n - 2]   superdiagonal of the matrix A
  306.  
  307.      b[0 .. n - 1]   right hand side, replaced by the solution vector x */
  308.  
  309.   size_t k;
  310.  
  311.   c[0] = d[0];
  312.  
  313.   if (n == 0)
  314.     {
  315.       return GSL_SUCCESS;
  316.     }
  317.  
  318.   if (n == 1)
  319.     {
  320.       b[0] = b[0] / d[0] ;
  321.       return GSL_SUCCESS;
  322.     }
  323.  
  324.   d[0] = e[0];
  325.   e[0] = 0;
  326.   e[n - 1] = 0;
  327.  
  328.   for (k = 0; k < n - 1; k++)
  329.     {
  330.       size_t k1 = k + 1;
  331.  
  332.       if (fabs (c[k1]) >= fabs (c[k]))
  333.     {
  334.       {
  335.         double t = c[k1];
  336.         c[k1] = c[k];
  337.         c[k] = t;
  338.       };
  339.       {
  340.         double t = d[k1];
  341.         d[k1] = d[k];
  342.         d[k] = t;
  343.       };
  344.       {
  345.         double t = e[k1];
  346.         e[k1] = e[k];
  347.         e[k] = t;
  348.       };
  349.       {
  350.         double t = b[k1];
  351.         b[k1] = b[k];
  352.         b[k] = t;
  353.       };
  354.     }
  355.  
  356.       if (c[k] == 0)
  357.     {
  358.       return GSL_FAILURE ;
  359.     }
  360.  
  361.       {
  362.     double t = -c[k1] / c[k];
  363.  
  364.     c[k1] = d[k1] + t * d[k];
  365.     d[k1] = e[k1] + t * e[k];
  366.     e[k1] = 0;
  367.     b[k1] = b[k1] + t * b[k];
  368.       }
  369.  
  370.     }
  371.  
  372.   if (c[n - 1] == 0)
  373.     {
  374.       return GSL_FAILURE;
  375.     }
  376.  
  377.  
  378.   b[n - 1] = b[n - 1] / c[n - 1];
  379.  
  380.   b[n - 2] = (b[n - 2] - d[n - 2] * b[n - 1]) / c[n - 2];
  381.  
  382.   for (k = n ; k > 2; k--)
  383.     {
  384.       size_t kb = k - 3;
  385.       b[kb] = (b[kb] - d[kb] * b[kb + 1] - e[kb] * b[kb + 2]) / c[kb];
  386.     }
  387.  
  388.   return GSL_SUCCESS;
  389. }
  390.